Valence Force Model for Phonons in Graphene and Carbon Nanotubes 
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Many calculations require a simple classical model for the interactions between sp^-bonded carbon 
atoms, as in graphene or carbon nanotubes. Here we present a new valence force model to describe 
these interactions. The calculated phonon spectrum of graphene and the nanotube breathing-mode 
energy agree well with experimental measurements and with ah initio calculations. The model does 
not assume an underlying lattice, so it can also be directly applied to distorted structures. The 
characteristics and limitations of the model are discussed. 



Graphene and carbon nanotubes are remarkable ma- 
terials, notable for both their fascinating properties and 
their technological promise In both contexts, it 

is often necessary to calculate the phonons for prob- 
lems where the use of ah initio methods is not feasi- 
ble. For graphitic systems, this has usually been ap- 
proached by approximating the force-constant matrix 
with terms coupling atoms up to some maximum distance 
0, [3, 0, H, @, 0, Si ■ This approach has many appealing 
features, but it has two important limitations. First, the 
terms in the force-constant matrix decay smoothly with 
distance between atoms so in practice it is necessary 
to truncate the expansion long before it has converged. 
Second, this approach is generally restricted to describ- 
ing phonons in the ideal crystal. It has required some 
ingenuity and inconvenience even to extend these mod- 
els to nanotubes, based on an idealized curved-graphene 
structure [1, i, S i] . 

In order to treat phonons in large low-symmetry sys- 
tems, such as rumpled graphene or bent nanotubes, one 
would like a model that explicitly gives the energy as 
a function of atomic positions, without reference to any 
underlying crystal structure. In principle one could use 
the general-purpose empirical interatomic potentials that 
are available for carbon, such as Ref. [T^]. But phonon 
applications typically require higher accuracy than such 
general-purpose models can provide. 

For diamond- and zincblende-structure semiconduc- 
tors, the problem was largely solved by the use of "va- 
lence force" models. These models use a smaller number 
of more complex terms, which may be more or less phys- 
ically motivated However, to date only one valence 
force model has been proposed for graphene [I^, ; and 
it explicitly references the graphene lattice, hindering ap- 
plication to distorted structures |14|. 

Here we present a new valence force model for sp^- 
bonded carbon. The model explicitly gives energy as a 
function of atomic positions, without reference to any un- 
derlying crystal structure. The only restriction is that the 
local geometry be consistent with sp^ bonding, i.e. three 
neighbors not too far from 120° degrees apart. Thus it 
can be directly applied not only to graphene, but to nan- 
otubes and fuUerenes, in relaxed or distorted configura- 
tions. We have tested the model for phonons in graphene 



and carbon nanotubes. We first describe the model itself, 
and the fitting procedure. We then present and discuss 
the phonon spectrum which is obtained after fitting the 
model parameters to selected experimental and theoret- 
ical data. Finally we discuss the overall accuracy and 
limitations, along with some related issues such as an- 
harmonicity. 

We write the energy as 

i,j<k£i 



3% ■ Vik X Vii 



(1) 



where 



being the atomic position vector 



of atom i, and the bondlength is rij = \vij\. In the 
summations, j € i means j runs over three neighbors of 
atom i, j < k € i means j and k are both neighbors of i 
(ordered to avoid double counting) , restriction j < k < I 
leaves only one possibility for the three neighbors of i, 
while restriction j ^ k < I gives three terms for each i. 
The bond length in the ground state of graphene is 



ro=0.142 nm; Sri 



SCi,jk — 



Vii 

3^ 



T'ij — fa ■ Wc further define 

TijTik 

X Vik + Vik X Vii + Vii X Vij 



niTik + TikTil + ruTi 



(2) 



where j, k, and I are the three neighbors of i. 

The first two terms in Eq. (U) represent the bond 
stretching stiffness B^^ and bending stiffness /3c, as in 
the Keating model [l5|. However, the form here avoids 
the large and unphysical anharmonicities of the Keating 
model. The third term f3y provides stiffness against out- 
of-plain vibrations. The fourth term Pre is motivated by 
bond-order potentials The fifth term (3p gives stiff- 
ness against misalignment of neighboring tt orbital. The 
last term Pre couples bond stretching and bond bending. 



2 



In fitting such a model, one typical chooses a set of data 
that it is desired to reproduce, and defines a weighted 
error which is to be minimized. As a straightforward 
test of the model and its ability to reproduce realistic 
phonon dispersions, we first try fitting to published LDA 
calculations 0, E^- The result is shown in Fig. [TJ [We 
follow the spectroscopic convention of reporting phonon 
energies in cm~^, where 1 cm~^ means hc/{l cm) « 0.124 
meV.] The rms error is only 22.6 cm~^, substantially less 
than the best previous fit to GGA dispersions using a 
valence force model with five parameters [isj . 




sponding phonon dispersion is shown as a solid curve in 
Figure [H The corresponding elastic constants are given 
in Table HIl (We equate in-plane elastic properties of 
graphene and graphite using the experimental layer spac- 
ing c = 6.7 A and volume per atom Vb = 3>/3roc/8.) 

Overall we consider the agreement in Figure [5] and Ta- 
ble [IT] to be quite good. The quality of the fit is a highly 
nonlinear function of the parameters, so there may be 
entirely different sets of parameters that give a similar or 
even better agreement with the same data. 



16001 



LU 




FIG. 1: Phonon band dispersions. Green dotted curves are 
LDA calculations 0]. Red solid curves are results of the model 
Eq. (m fitted to these LDA calculations. The corresponding 
parameters are given in [l^ . 



By giving more weight to one feature or another in 
the fitting, it is straightforward to improve the descrip- 
tion of e.g. the acoustic branches at the cost of worsen- 
ing the optical branches. However, regardless of how we 
weighted the data, we could not reproduce the dips in 
the highest phonon bands at F and K while keeping a 
reasonable overall dispersion. This issue was also men- 
tioned in [isj . Electron-phonon interactions are known 
to affect phonon dispersions even in bulk semiconduc- 
tors ; and such interactions are particularly important 
for the highest bands of graphene near F and K due to 
the Kohn anomaly Thus we cannot expect to de- 

scribe these dips with short-range classical interactions. 
It would therefore seem logical to fit the bands away from 
F and K, and accept that the model gives energies that 
are too high for the top bands at those points. However, 
because the optical phonon energy at F is a widely used 
reference, we have chosen to fit this point accurately. 

We find that the Poisson ratio calculated with our 
model fitted to the LDA calculations alone \s v = 0.4, 
much less than the experimental value of = 0.17. This 
suggests that elastic properties should be included in the 
fitting. Also, the experimental and theoretical data are 
not in perfect accord. We have therefore chosen to fit a 
mixture of published experimental phonon data, ah initio 
phonon calculations, and elastic constants. The result- 
ing parameter values are listed in Table [H and the corre- 



FIG. 2: Phonon band dispersion, comparing fitted model with 
experimental data and ab initio calculations. Red solid curve 
is our model, Eq. ([T}, with the parameters given in Table 
m Green dotted curve is an LDA calculation [3. Blue sym- 
bols are experimental data: electron energy loss spectroscopy 
(EELS) from Refs. [3, and [23| (respectively squares, 
diamonds, and filled circles), neutron scattering from Ref. [2l| 
(open circles), and x-ray scattering from Ref. [221 (triangles) . 
Data for Refs. O, [M S EH, HI are taken from Ref. 



TABLE I: Parameters of the model Eq. H} used in Fig. [2] 
based on best fit to the experimental data and LDA calcula- 
tions. Units are eV. 





/3c 




/3.2 






18.52 


4.087 


1.313 


4.004 


0.008051 


4.581 



The longitudinal and transverse sound velocities {v 
doj/dq at q = 0) within our model are 



81 



4 12/3^1 4- 27(3c - 6/3^2 - 18/?^ 



Mcvl = McV^ + - {Prl + Pr2) 



(3) 



where Mc is the mass of a carbon atom. The model 
velocities vta ~ 13.3 km/s and vla ~ 21.2 km/s are very 
close to the experimental values of vta ~ 14 km/s and 
WLA ~ 24 km/s [l3|- The elastic constants are related to 
the sound velocities: VqCqq = AIcv'^ and VqCu — Mcv\ 
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It is often convenient to have analytic expressions for 
the phonon energies at symmetry points (e.g. for verify- 
ing a numerical implementation). From Eq. ([1]), 



(4) 



where the index i runs over (rl, c, v, r2,p, re) and coeffi- 
cients ai are given in Table IIIII 

Turning from graphene to carbon nanotubes, we cal- 
culate the radial breathing mode (RBM) for tubes of dif- 
ferent diameter and chirality. This mode corresponds to 
a radial stretching or compression of the tube. The mode 
emerges from the lowest-energy acoustic phonon modes 
in graphene. The RBM acquires finite energy at zero 
wavcvector due to the nanotube curvature, with a simple 
K, \/d scaling of energy with diameter. As a result, RBM 
measurements are widely used to identify the diameters 
of single- walled carbon nanotubes [2^. 

For a given tube, we first relax the atomic positions 
and allow the lattice constant to adjust to minimize the 
energy. We then calculate the RBM energy. The results 
for all tubes in the diameter range from 0.5 to 4.0 nm are 
shown in Fig.[3l Simple scaling arguments based on con- 
tinuum elasticity suggest that RBM energies should scale 
with diameter as Hujrbm — A/d. Experimental data 
arc typically fitted with the phenomcnologically adapted 
form hionBM = A/d + B. For tubes of d=l nm, exper- 
imental phonon energies A-l-B ar e repor ted in the range 



226-248 cm-^ [27|, |28|, |2g|, 

calculations suggest A-\-B 



31, 32, 3 



while ab-initio 

■ 1 



=234 or 226 cm-^ (271,133. The 

constant off-set was reported in the range from B = —6 
cm~^ to B = 27 cm"^. Recently, it was reported [s^ that 
a non-zero offset B is caused by the interaction with a 
substrate, while for freely suspended nanotubes B should 
be zero. 

Within our model, the RBM mode shows accurate A/d 
scaling independent of chirality, with B k, Q and A « 225 
cm~^ (where d is in nm) as shown in Fig. [3) From the 
theory of elasticity, A = 2hvL, which gives A k, 225 
cm^^ for the parameters of Table HI in accord with the 
numerical result. The model is in good agreement with 
the most recent experimental [s^ and theoretical (3^ 
values of A =-221 and A + B =226 cm^^ respectively. 

In general, a valence force model will have some an- 
harmonicity. Since wc have not attempted to fit exper- 
imental or ab initio anharmonicities, any anharmonic- 



TABLE II: Elastic constants from the experiment and the 
model (in GPa). Note a relation among elastic constants 
[3, [2^ for hexagonal symmetry: Ces = (Cn — Ci2)/2, 

/^ = Cl2/Cll. 





Cii 


C12 


Gee 


V 


experiment [24] 


1060 


180 


440 


0.17 


model 


1024 


210 


407 


0.20 



ities are likely to be unphysical. It is therefore desir- 
able to minimize the anharmonicity in the model, and 
the form of Eq. ([1]) is designed with this in mind. One 
measure of anharmonicity is the Gruneisen parameter 
7 = —(2u))~'^{dLLi/de), which represents the fractional 
shift in phonon frequency uj when the crystal is subjected 
to a strain e in all directions. For the doubly degen- 
erate phonon mode in graphene, our model gives 
^E'lg ~ —0.2. This is much smaller in mag nitude than 
the experimental value of ^E^g ~ 2.0 13l. l34|. confirming 



that our model is relatively harmonic in this respect. 

For nanotubes, we have another form of anharmonicity, 
the phonon shifts due to bending of the graphene sheet. 
We have calculated the shifts in LO and TO phonons 
relative to graphene. The TO mode shift is less than 
12 cm~^/d^ in our model (where d is in nm), and the 
LO mode and the LO-TO splitting are even less. Ex- 
perimental shifts are four times larger in semiconducting 
nanotubes (ssj . confirming that our model successfully 
minimizes any unintended anharmonicities. 

In conclusion, we have developed a valence force model 
applicable for sp^-carbon based structures. Our model 
gives a good fit of the graphene phonon dispersion and 
elastic constants, and describes well the RBM energy of 
nanotubes. The model also avoids the unphysical strong 
anharmonicities that occur in some valence force models. 
Most importantly, in contrast to other phonon models 
for sp^-bonded carbon, Eq. ([T]) makes no reference to 
an underlying lattice, so it can be directly applied to 
distorted geometries. 

We gratefully acknowledge N. Marzari, O. Dubay, and 
D. Kresse for providing data for Figs. [2] and (U and 
W. A. Harrison and A. Jorio for helpful discussions. 
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FIG. 3: Radial Breathing Mode (RBM) energy as a function 
of tube diameter (red open circles) along with the best fit 
huRBM = 224.6 cm~^/d (black solid curve). The inset shows 
the same results versus inverse diameter. 
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